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Abstract 

The formation of a single bubble from an orihce in a solid surface, submerged in an in¬ 
compressible, viscous Newtonian liquid, is simulated. The hnite element method is used to 
capture the multiscale physics associated with the problem and to track the evolution of 
the free surface explicitly. The results are compared to a recent experimental analysis and 
then used to obtain the global characteristics of the process, the formation time and volume 
of the bubble, for a range of orihce radii; Ohnesorge numbers, which combine the material 
parameters of the liquid; and volumetric gas how rates. These benchmark calculations, for 
the parameter space of interest, are then utilised to validate a selection of scaling laws found 
in the literature for two regimes of bubble formation, the regimes of low and high gas how 
rates. 
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1. Introduction 


The controlled production of small gas bubbles is of critical importance to many processes 
found in the chemical, petrochemical, nuclear, metallurgical and biomedical industries. Sub¬ 


sequently, a vast body of research exists on the subject (Kumar and Kuloor, 1970 Clift et 


al., 1978 

Kulkarni and Joshi 

2005 

Yang et al., 2007 

). Due to the complex mechanisms 

involved in this fundamental process 
erature has focused on the case of f 
formation site, namely a submerged 

, the computational, experimental and theoretical lit- 
^enerating a bubble by pumping gas through a single 

vertical nozzle (Longuet-Higgins et ah, 1 

991 Oguz and 

Prosperetti, 1993; Wong et al. 

1998 

Thoroddsen et al. 

2007) or, as is the case in this work. 

an upward facing orihce in a submerged solid surface (5 

jhang and Shoji| 2001 

i Badam et al. 

2007 Gerlach et ah, 2007 Das et al. 

2011). 


Some authors have considered the case of inhating a bubble under a constant gas pressure 
with the volumetric gas how rate varying as the bubble inhates (fDavidson and Schuler 


Satyanarayan et ah, 1969 La Nauze and Harris, 1972). Others have considered 
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the case where the volumetric gas flow rate is determined by the difference between some 
ambient pressure away from the bubble and the gas pressure in a chamber that is connected 


to the formation site (Khurana and Kumar, 1969 Oguz and Prosperetti, 1993). However, 


the most popular method, and the method employed here, is to apply a constant volumetric 


gas flow rate through the formation site and allow the gas pressure to vary in time (Wong 


et ah, 1998 Jamialahmadi et ah, 2001; Corchero et ah, 2006 Gerlach et ah, 2007). 


1.1. Experimental Observations 


The majority of experimental studies examine a continuous chain of bubbles (Zhang and 


Shoji, 2001 Badam et al.| 2007 Das et ah, 2011| . A bubble grows whilst attached to the 
formation site as gas is pumped through the site at a constant flow rate. As the volume 
of the bubble increases, the influence of buoyancy becomes more important. The bubble 
seeks to minimise its surface area at a given volume and so a ‘neck’ develops in the bubble 
as the longitudinal curvature of the free surface changes sign at some point just above the 
three phase solid-liquid-gas contact line. The difference in pressure between the base and 
the apex of the bubble then drives the thinning of this neck as the free surface begins to 
pinch-off (see Figure [^. This leads to the eventual break up of the bubble into two parts; 
a new bubble is released and rises away from the formation site under buoyancy, whilst the 
residual bubble, that is left attached to the formation site, begins to grow. This process 
repeats itself, thus producing a chain of bubbles. 

The aim of many of these studies has been to And the global characteristics of the 
formation process, such as the frequency of formation and the volumes of the bubbles that 
are formed, for a given set of material parameters (e.g. liquid density and viscosity), design 
parameters (e.g. orifice radius, wettability of substrate) and regime parameters (e.g. gas 
flow rate). These studies provide potentially useful material for comparing theoretical results 
with experiments. However, global characteristics of bubble formation in a chain accumulate 
several complex phenomena, and it is desirable to have a more detailed picture of each of 
them and of their interaction. 

Theoretical and early computational studies focused on the behaviour of a single bubble 
which grows from some initial state up until the break up of the free surface is approached 


'Longuet-Higgins et ah 

1991 

Oguz and Prosperetti 

1993 

Wong et al. 

1998 

Xiao and Tan 


2005). The mathematical difficulties of handling the topological change associated with the 


complete break up of the bubble prohibited any further progress and so, in this case, the 
formation time period and the volume of the bubble above the point of minimum neck radius 
are of interest and have been investigated. Here, the break up is assumed to be a local effect 
that does not affect the global dynamics of the formation process. 

Experimental and theoretical studies identified three regimes for the formation of bubbles 


under a constant gas flow rate (McCann and Prince, 1971 Zhang and Shoji, 2001 Badam 


et al.| 2007). For a given set of material and design parameters, the first of these regimes 
occurs at small gas flow rates and is known as the ‘static’ regime. In this regime, the volume 
of the bubbles formed is independent of flow rate, and therefore a decrease in the flow rate 
results in an increase in the formation time. Consequently, it is not possible to produce a 


bubble with a volume smaller than this limiting volume. Fritz (1935) suggested that this 
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limiting final volume is proportional to the orifice radius. This, as well as other scaling laws 
are considered in section |5l 

At greater flow rates, bubble formation enters the ‘dynamic’ regime. It is now the 
formation time that approaches a limiting value, resulting in an increase of bubble volume 
Vd with flow rate Q. In this regime, some authors used spherical bubble models to propose 


scaling laws for the bubble volume. Davidson and Schuler (1960b) and Oguz and Pr osperetti 


(1993) proposed that oc Qs for bubble formation in an inviscid liquid, whilst, Wong et 
.. 


(1998) suggested that Vd oc Qj in the case of a highly viscous liquid. 

'i^inally, under even greater flow rates, bubble formation enters the ‘turbulent’ regime. 
Here the motion is chaotic with successive bubbles coalescing with each other above the 
formation site. The transition between these three regimes is dependent on all of the material 
and design properties involved. For instance, in the case of an orihce, where the liquid phase 
only partially wets the solid surface and so the contact line is free to move along the solid 


surface, the wettability of the substrate is important to the formation process (Lin et ah 


1994 Byakova et ah, 2003 Gnyloskurenko et ah, 2003 Corchero et ah, 2006). 


1.2. Theoretical Progress 

Under sufficiently small gas flow rates, the initial growth of a single bubble may be 
accurately modelled as a quasi-static process, an approach that has also been utilised for 
the description of the formation of liquid drops (Fordham, 1948; Thoroddsen et al., 2005). By 


assuming that the fluid velocities are negligible, one can derive the Young-Laplace equation, 
which balances the difference in the gas and hydrostatic pressures with the capillary pressure. 
This can be solved to hnd the free surface prohle of a bubble at a particular volume. By 
stipulating a small increase in the bubble volume, a series of successive prohles can be found 
which are in very good agreement with experiments that describe the initial evolution of a 


bubble (Marmur and Rubinf 1 

973[ Longuet-Higgins et al 

1991| Gerlach et al. 

2005| |Lee 

and Tienl,|2009t|Vafaei and We 

n, 2010t Vafaei et ah, 2010, 

font Lee and Yang, 2 

3121 Lesage 

and Marois, 2013). 


However, once the neck forms and the pinch-off process begins, the liquid velocities 
adjacent to the point of minimum neck radius are no longer negligible. Due to the dynamics 
associated with the relatively large liquid velocities involved at this stage, the quasi-static 
approach is no longer valid and thus, in general, one cannot expect to accurately predict the 
hnal volume of a newly formed bubble this way. 

In an attempt to describe some aspects of the dynamic problem, early mathematical 
models of single bubble formation under a constant gas flow rate were based upon global 


force balances. The hrst ‘one-stage’ models for highly viscous (Davidson and Schuler, 1960a) 


and inviscid liquids (Davidson and Schuler, 1960b), which involved the bubble growing 


spherically, were developed into various ‘two-stage’ models for highly viscous (Ramakrishnan 


et al., 1969 Gaddis and Vogelpohl, 1986| and inviscid liquids ([Wraith 1971; Buyevich and 


Webbon, 1996) by adding a detachment stage to the spherical expansion stage. Once the 


spherical bubble had reached a certain volume, it translates away from the formation site 
but remains in contact with it via a cylindrical column of gas. The bubble is said to detach 
when the column reaches a certain height. In order to incorporate intermediate values 
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of viscosity, some authors solved a modified Rayleigh-Plesset equation, which describes a 
spherical bubble oscillating under pressure, by assuming that each element on the free surface 
is the corresponding sphere with the same longitudinal radius of curvature. From this, the 
mean gas pressure could be found that was then used to solve an equation for the upwards 


translation of the bubble (Pinczewski, 1981 Terasaka and Tsuge, 1990, 1993). Although 


these semi-empirical models are seen to give good agreement in certain regimes (Kulkarni 


and Joshi, 2005), to accurately describe the whole parameter space of interest, the use of 


complex computational techniques is needed. 

1.3. Computational Approaches 

To solve the unsteady free-boundary problem of bubble formation from an orihce or 
nozzle, subject to the forces of gravity, inertia, viscosity and capillarity, the development of 
computational fluid dynamics techniques is required. Much of the early work in this direction 
was concerned with the axisymmetric generation of a single bubble in either the inviscid 


'Marmur and Rubin 

1976 

Hooper 


Tan and Harris, 


Oguz and Prosperetti 


1993 Xiao and Tan, 2005; Higuera and Medina, 2006||) or highly viscous flow regime (Wong 


et al., 1998 Higuera, 2005), where the boundary integral method can be used to reduce the 


problem’s dimensionality to one. 

Since then, the majority of work on bubble generation has focused on the influence of the 
material, design and regime parameters on the global characteristics of the flow using the 


volume-of-fluid method 

(Ma et al. 

2012 

), and its improved variants which utilise level-set 

methods ( 

Buwa et al. 

, 2 

007 

; ( 

Jerlach et al., 2007 

; 

Jhakraborty et al., 

2009 

Ohta et al. 

2011 

Chakraborty et al. 

2011 

Albadawi et al. 

201f 

1). Due to the simple manner in which 


topological changes are ‘automatically’ handled, such techniques have proved successful at 
describing many features of the bubble formation phenomenon, including the wake effect 
of a preceding bubble on subsequent bubbles in a chain, where the size of forming bubbles 
are assumed sufficiently large and hence the details of how the topological transition (i.e. 
the break up of the bubble) takes place are relatively unimportant. However, as noted in 
Eggers and Villermaux ( 2008[ ), there is no guarantee that the physics associated with this 
transition has been properly accounted for. This aspect becomes more important when the 
size of forming bubbles is small. 

Despite the advantages of volume-of-fluid based methods, as a relatively simple way of 
handling the global dynamics, it is well known that this class of numerical techniques are not 
well suited to resolving the multiscale physics which becomes critical in ‘singular’ flows, i.e. 


those in which liquid bodies coalesce (Hopper, 1990; Eggers et al., 1999) or divide (Rayleigh 


1892 Eggers, 1997). For instance, for a millimetre-sized bubble, experiments are able to 


resolve the minimum neck radius down to tens of microns (Burton et al., 2005 Thoroddsen 


et al., 2007; Bolanos-Jimenez et ah, 2008, 2009), whilst numerical methods have thus far 


failed to capture these scales and often artihcially truncate the simulation far above the 
scales which are still well within the realm of continuum mechanics. These problems have 


been highlighted in the work of Hysing et al. (2009) 


These dehciencies can be addressed by using the hnite element method which has pre¬ 


viously been used to capture inherently multiscale flows in dynamic wetting (Wilson et al. 
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2006 Sprittles and Shikhmurzaev, 2013) and in the coalescence of liqnid drops (Sprittles 


and Shikhmurzaev, 2012c). Moreover, it can easily be applied to the bubble formation 


phenomenon so that (a) the parameter-space of interest to technologically-relevant pro¬ 
cesses can be investigated, safe in the knowledge that all scales of the problem have been 
accurately resolved, and global characteristics of the flow, such as the formation time and 
bubble volume, can be extracted and (b) the actual pinch-off of the bubble can be studied 


and compared to various theories proposed in the literature for ‘singular’ flows (Burton et 

al.l 2005 Gordillo et ah 

2005 Gordillo and Perez-Saborid, 2006 |Thoroddsen et ah, 2007 

Bolanos-Jimenez et ah 

2008t Quan and Hua 

2008 IGordillol 2008; Bolanos-Jimenez et ah 

2009[ Gekle et ah, 2009 

Fontelos et al.| 2011 

). Historically, aspects (a) and (b) have been 


considered by separate scientihc communities, but in this work, techniques are used which 
are able to resolve both features. Here, the focus will be on aspect (a), i.e. a parametric 
study, with a forthcoming publication studying aspect (b) in detail. 

The outline of the paper is as follows. In section the problem formulation is given 
and three dimensionless parameters are identified. Section then contains the range of 
physically-realistic values of these parameters that are to be used to investigated each 
parameter’s influence on the bubble formation problem. A brief description of the finite 
element computational framework used is given in section The results are presented in 
section where a comparison with an experiment can be found together with the bench¬ 
mark calculations, which describe the influence of the dimensionless parameters on the global 
characteristics of the bubble formation process for the parameter space of interest. These 
calculations are then used to validate various scaling laws found in the literature. Finally, 
some concluding remarks can be found in section 


2. Problem Formulation 

Consider a smooth, horizontal, stationary, impermeable solid surface submerged in a 
quiescent, incompressible, viscous Newtonian liquid of constant density p and dynamic vis¬ 
cosity p. The solid surface has a circular orifice of dimensionless radius through which 
an inviscid gas is pumped at a constant dimensionless volumetric flow rate Q to form a 
single bubble. The characteristic velocities in the gas and the size of the resulting bubble 
are assumed to be sufficiently small so that any spatial non-uniformity of the gas pressure 
in the bubble can be neglected. 

The following scales are dehned for length L = \/a/pg, velocity U = a/p, pressure 
cr/L = ^/pag, time L/U = p/y/f^ and flow rate L‘^U = a‘^/ppg, where a is the surface 
tension of the gas-liquid interface and g is the acceleration due to gravity. In other words, the 
scales are based on a regime in which viscous, capillarity and buoyancy forces are of similar 
magnitude, so that the Bond number {Bo = pgL'^/a) and capillary number {Ca = pU/a) 
are unity. From here on, all quantities are dimensionless unless stated and where dimensional 
quantities do appear, they are denoted with bars. 

The incompressible flow of the liquid is governed by the dimensionless Navier-Stokes 
equations, 

V • u = 0, (la) 
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(lb) 


—+ (u-V)u = 0/i2 (V-P-e,), 

where u = u(r, t) and p = p(r, t) are the liquid velocity and pressure, r is the spatial 
coordinate vector, t is time, P = —pi + Vu + (Vu) is the stress tensor, I is the metric 
tensor and is a unit vector in the opposite direction to gravity. The Ohnesorge number 
is Oh = fi/y/paL which, by substituting the length scale, is. 


Oh 


h9 


1/4 




The kinematic boundary condition on the free surface, that states the fluid particles that 
form the free surface remain there for all time, is. 


dl 

dt 


+ U-V/ 


0 , 


( 2 ) 


where /(r, t) = 0 dehnes the a-priori unknown free surface. The balance of forces acting on 
an element of the free surface from the liquid and gas phases and neighbouring surface ele¬ 
ments manifests itself through two dynamic boundary conditions, the normal and tangential 
stress conditions, 

Pg n ■ P ■ n = V ■ n, (3) 

n ■ P ■ (I — nn) = 0, (4) 

respectively, where Pg = Pg{t) is the spatially uniform gas pressure, which is determined as 
part of the solution, and n is the unit normal to the free surface that points into the liquid 
phase (see Figure [^. The combined no-slip and impermeability condition is applied on the 
solid surface, 

u = 0. (5) 


The fluid flow is assumed to be axisymmetric about the vertical axis and therefore the 
problem is considered in the (r, 2 ;)-plane where r and 2 ; are the respective radial and vertical 
coordinates of a cylindrical coordinate system (see Figure [^. The velocity and pressure 
helds of the liquid can therefore be rewritten as u{r,z,t) = u{r,z,t)er + w{r,z,t)ez and 
p(r, z, t), respectively, where is the unit coordinate vector in the r direction. The a-priori 
unknown free surface is now dehned as f{r,z,t) = 0. 

The conditions of axial symmetry and smoothness of the velocity held at the axis of 
symmetry are given by. 


u = 0, 



( 6 ) 


whilst in the far held, there is. 


■u, w —)■ 0, 


r'^ + z‘^ ^ 00. 


( 7 ) 


The normal stress boundary condition ^ on the free surface is itself a diherential equa¬ 
tion determining the free surface shape and as such requires its own boundary conditions. 
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Figure 1: A sketch of the flow domain in the (r, ; 2 )-plane. 


Firstly, it is assumed that the liquid phase perfectly wets the solid surface and so the three 
phase solid-liquid-gas contact line remains pinned to the rim of the orifice for all time, 

f{rc,0,t) = 0 , (8) 

where the contact angle 9 = 9{t) at the pinned contact line is then determined as a part of 
the solution. Secondly, the free surface is smooth at the bubble apex {r = Q, z = H) and so, 

V/ ■ e, = 0, (9) 

is applied there, where H = H{t) is the height of the bubble to also be determined as part 
of the solution. 

For the initial-value problem, the initial velocity held and the initial shape of the bubble 
must be prescribed. It is assumed that the liquid is initially at rest, 

u{r, z, 0) = 0, w{r, z, 0) = 0, (10) 

and that the bubble is initially a spherical cap with volume Vi. The initial height of the 
bubble Hi = H{0) is the real solution of nHi (3r^ + Hf) = 6Vi and the radius of the sphere 
is dehned by Ri = {Hf + r^) /2Hi. The initial free surface shape is then given by, 

+ {z + Ri — Hif = 0 < 2 : < Hi. (11) 

In other words, by dehning the radius of the orihce and the initial volume of the bubble, the 
initial shape is fully specihed. 
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To conclude the problem formulation, an equation governing the volume of the bubble is 
required in order to employ a constant volumetric gas flow rate. The volume of the bubble 
is given by, 

V{t) = Vi + Qt, (12) 

and so if the free surface is described by the function r = h{z, t) where /(r, z, t) = h{z, t) — r 
then, 

nz=H 

V(t)=7i / h‘^{z,t)dz. (13) 

Jz=0 

The topological change associated with the break up of the bubble, as the simply con¬ 
nected gas phase becomes multiply connected, has not been accounted for in this problem 
formulation and therefore any simulations will stop short of the break up of the bubble. 
Each simulation runs until the minimum neck radius = TtoZ, where the tolerance rtoi is 
discussed in section HI 


3. Parameter Regime 


The problem can be characterised by the three dimensionless parameters identified in 
the problem formulation above; the orifice radius r^, the Ohnesorge number Oh and the 
volumetric gas flow rate Q. 

To estimate the parameter regime of interest, consider three typical Newtonian liquids; 
water (/i = 1 mPa s, p = 998.2 kg m“^, a = 73 mN m“^), a silicone oil {p = 0.01 Pa s, 
p = 800 kg m“^, cr = 20 mN m“^) and glycerol {p = 1.4 Pa s, p = 1200 kg m“^, a = 
60 mN m“^). The respective length scales and Ohnesorge numbers for these three liquids 
are then L = 2.73 mm, 1.60 mm and 2.26 mm and Oh = 2.24 x 10“^, 6.26 x 10“^ and 3.49. 

In order to investigate the full effects of these parameters, the system shall be examined 
using the orifice radii Vc = 0.1 and 1, since it is the generation of bubbles from orifice 
radii of the order of millimeters and below that are of interest; the Ohnesorge numbers 
Oh = 2.24 X 10“^, so that the case of water can be examined explicitly, and also Oh = 10“^, 
10“^, 1 and 10, to cover the full range of liquid viscosity; and flow rates 10“® < Q < 0.5 
for Tc = 0.1 and 10“^ < Q < 15 for Vc = 1. Each simulation will be characterised by the 
notation (uc. Oh, Q). 

The initial volume of the bubble Vi is chosen to be sufficiently small enough so that, 
even though a smaller initial volume would slightly increase the formation time t^, it would 
not alter the volume of the newly formed bubble Vd- It was seen for all simulations that 
an initial contact angle of 0(0) = 37r/4 is suitable to fulhll this criterion and so the initial 
volumes for the cases of rc = 0.1 and 1 are Vi = 6.89 x 10“^ and 6.89 x 10“^, respectively. 

The length and velocity scales used here are different from those which the reader may 
have encountered in similar works. A detailed guide on how to rescale the results presented 
in this manuscript can be found in Appendix A[ 
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4. A Computational Framework for Bubble Formation 


Due to the mathematical complexity of the bubble formation problem, where gravita¬ 
tional, viscous, inertial and capillarity forces are all present in an unsteady free-boundary 
problem, a computational method is required to solve the dimensionless system of equations 
0-1111. The hnite element computational framework used here is based on a formula¬ 
tion described in detail in |Sprittles and Shikhmurzaev (2012a, 2013), where a user-friendly 
step-by-step guide to its implementation can be found alongside numerous benchmark cal¬ 
culations. This method has already been used to simulate both drop impact phenomena 


’Sprittles and Shikhmurzaev 

2012b 

) and, more recently, two-phase coalescence processes 

’ Sprittles and Shikhmurzaev, 

2014 

)• 

Consequently, the method is only brieffy described here 


and the reader is referred to the aforementioned references for further details. 

The problem is formulated in an infinite domain (r > 0, z > 0) with boundary conditions 
([^ prescribed in the ‘far-field’. In order to apply the finite element method, the fiow domain 
must be truncated to a finite size. Therefore, a side boundary {r = Q < z < Z) and a 
top boundary {Q < r < z = Z) are constructed on which amended ‘passive’ boundary 
conditions are applied. To ensure these amended conditions are indeed passive, R and Z 
are set large enough so that any increase in their values does not affect the growth of the 
bubble. This is achieved when R is at least five times the maximum width of the bubble 
and Z is at least four times the maximum height of the bubble. 


As is standard in the finite element method (Gresho and Sani, 1999), the system of 


equations is discretised over a set of nodes that are distributed throughout the fiow domain, 
whilst a finite number of non-overlapping triangular elements are tessellated about these 
nodes to form a computational mesh (see Figure [^. The V6-P3 elements are used to 
approximate velocity quadratically and pressure linearly thus satisfying the Ladyzhenskaya- 
Babuska-Brezzi condition (Babuska and Aziz, 1972[ ), whilst curved elements allow the free 
surface to also be captured with a quadratic approximation. 

An arbitrary Lagrangian-Eulerian mesh design based on the ‘method of spines’ ([Kistler 


and Scriven, 1983) is used where nodes are distributed along the free surface to track its 


evolution explicitly, whilst the nodes in the bulk move in a prescribed manner that ensures 
the elements remain non-degenerate. 

Figure shows the evolution of the computational mesh for the simulation shown in 
Figure Figure shows that the finite element method is ideally suited to deal with the 
bubble formation process as the mesh can be refined by placing smaller elements adjacent 
to the bubble where additional precision is required. The size of the elements then increase 
away from the bubble, where fiow variables vary on a larger spatial scale, so that the problem 
remains computationally tractable. 

To make clear the arrangement of these smaller elements adjacent to the free surface. 


Figures and show a close up of the bubble. To accurately capture the dynamics of 
bubble growth, i.e. the shape of the detaching bubble as well as the pinch-off of the neck, 
it is found that the mesh adjacent to the neck of the bubble requires more refinement than 
the mesh above the neck. 

All nodes move along straightened spines apart from those nodes in the pinch-off region. 
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(a) The entire flow domain corresponding to the 


computed prohle in Figure 3b 



(b) A close up of the computational 
mesh around the bubble from Figurej^ 



(c) A close up of the computational 
mesh adjacent to the bubble corre¬ 
sponding to the computed profile in 


Figure 3d 


Figure 2: The evolution of the computational mesh for the simulation shown in Figure 


Here, the spines are designed using a bipolar coordinate system with a focus at the contact 
line. This arrangement is ideal to capture the dynamics of the neck of the bubble as it thins. 
It was also designed with future research in mind, to capture the dynamics associated with 
a contact line that is free to move along the solid surface, i.e. the case of a liquid phase that 
only partially wets the solid surface. In addition, as will be seen later (see Figure [I^, under 
certain conditions, the liquid phase may enter the mouth of the orihce, and this design, in 
contrast to simpler ones, ensures the free surface can still be captured accurately when this 
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situation arises. 

Each one-dimensional element placed along the free surface consists of three nodes, with 
the two end nodes being shared with the two neighbouring surface elements. At each time 
step, the smaller surface elements in the pinch-off region and the larger surface elements 
around the rest of the bubble are equally spaced, respectively. For all simulations, at least 
125 nodes (62 surface elements) were used to track the free surface. A majority of these 
nodes, 81, were located in the pinch-off region. In some cases of very large flow rates, 
simulations used 153 nodes (76 surface elements) due to the increase in the size of the 
bubble and the necessity to include an additional region of the mesh, within the pinch-off 
region and adjacent to the contact line, to adequately capture the deformation of the surface 
that occurs there at the beginning of a simulation. Further increases in the number of nodes 
is seen to have a negligible impact on the global characteristics of the flow and, remeshing, 
which can have an adverse effect on accuracy, is not required due to the well-designed 
meshing strategy. 

A simulation runs until the minimum neck radius = 'r'toi, where the tolerance rtoi = 
5 X 10“^ Tc- Since it is the global characteristics of the bubble that are of interest here and, 
as stated previously, the problem should remain computationally tractable, more rehnement 
is not added and the simulation stops at this point. As = rtoi, the distance between 
adjacent nodes on the free surface surrounding the point of minimum neck radius for the 
largest bubble reported is 5 x 10“^ or 10“^ for the cases of Tc = 0.1 and Tc = 1, respectively. 

A second-order backward differentiation formula is used to calculate all derivatives with 
respect to time. A constant time step is used up until the bubble forms a neck and begins to 
pinch-off. At least 100 time steps were used before the neck is formed. This was seen to be 
small enough so that any decrease in the time step did not sufficiently affect the accuracy 
of the solution. As the neck thins, around 500-700 time steps are used as the time step 
decreases logarithmically with the minimum neck radius. As rtoi, the time step is 

O {IQ-^yO (10-5) for r, = 0.1 and O (lO'^) for r, = 1. 


Since the free surface of the bubble is the only boundary of the flow domain that is 
not hxed, as the bubble continues to grow, the size of the flow domain decreases. In order 
to ensure incompressibility in the flow domain, a ‘reference’ pressure of the liquid must be 
prescribed at a particular point for the entire simulation. This arbitrary reference pressure 
should not affect the growth of the bubble and so the straightforward choice would be to 
base this pressure on the assumption that the pressure held along the side boundary should 
remain hydrostatic, p{R, z, t) = —z, for the entire simulation. When p(i?, 0, f) = 0 is applied 
as the reference pressure, the point (i?, 0) behaves like a point sink which affects the velocity 
held close to the bubble. As this is far from acceptable, p(i?, Z, t) = —Z is applied instead. 
With the values of R and Z stated above, the pressure held along the side boundary remains 
hydrostatic for the entire solution. 

In order to validate the numerical platform, various simulations of a Rayleigh bubble 
were carried out and the results compared to those obtained from the Rayleigh-Plesset 
equation (Plesset and Prosperetti, 1977). The simulations were seen to converge to the 


analytic solution with increased spatial and temporal resolution, as expected. 
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fc) tp = 2.340 s and t = 2.927 s. 


(d) tp = 3.550 s and t = 4.137 s. 


Figure 3: The superposition of computational free surface profiles and experimental im¬ 
ages for the case {rc,Oh,Q) = (0.106,2.24 x 10“^, 5.11 x 10“®). The experimental images 


are reproduced from Di Bari and Robinson (2013). The dimensional experimental 4 and 


simulation t times are given with each subhgure. 
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5. Results 

The influence of the three dimensionless parameters, the orifice radius Tc, the Ohnesorge 
number Oh and the gas flow rate Q, on the global characteristics of bubble formation, in 
particular on the dimensionless formation time td and volume V^, can now be investigated. 
The formation time td is taken to be the time of the final solution, as ^ ^tou whilst the 
volume Vd of the bubble that is formed is the volume of the bubble above the point on the 
free surface which represents the minimum neck radius at the formation time. Where the 
free surface is displayed in a figure, the solid surface will also be shown as a horizontal black 
line at r > Tc, z = 0. 


Figures (3a)-(3d) comprise of an 


5.1. A Typical Case of Bubble Formation 

Figure shows a typical case of bubble formation, 
experimental image, reproduced from a recent paper of Di Bari and Robinson (2013), with 
the corresponding free surface profile, computed by the numerical platform, superimposed 
on top. Since Figure |^shows the entire cross-section of the bubble, the dimensional abscissa 
is given by x rather than f. The experimental images show the formation of an air bubble 
in water (p = 998.2 kg m“^, /i = 10“^ Pa s, cr = 73 mN m“^) by applying a dimensional 
volumetric gas flow rate oi Q = 2.78 mm^ s“^ through an orifice of dimensional radius 
fc = 0.29 mm in a submerged solid surface. There is very good agreement between the 
simulation and the experiment. 

With g = 9.81 m s“^ and using the scales identified in the problem formulation, the 
material parameters of the liquid give rise to the length scale L = 2.73 mm, velocity scale 
f/ = 73 m s“^, pressure scale cx/L = 26.7 Pa, time scale L/U = 3.74 x 10“® s and flow rate 


scale LfU = 5.44 x 10 ^ m^ s The experiment is then classified as the dimensionless case 
of {rc,Oh,Q) = (0.106,2.24 x 10-3,5.11 x lO"®). 

The dimensional experimental times tg that accompany each image are 9.60 x 10“^ s, 
1.380 s, 2.340 s and 3.550 s, respectively. However, the initial configuration at 4 = 0 is 
unclear as that result was not published and so, in order to compare the simulation to the 
experiment, the simulated free surface profile that best matched the fourth experimental 
image (see Figure 3d) was selected. This profile occurred at the dimensional simulation 
time f = 4.137 s, and so the criterion 7=4 + 0.587 was used to And the corresponding 
simulation time of the first three experimental images. The dimensional simulation times of 
the four computed solutions are therefore t = 0.683 s, 1.967 s, 2.927 s and 4.137 s. 

At 7 = 0 s, the simulated bubble starts as a spherical cap of volume 1.677 x 10“^ mm^. 
Initially, the bubble grows spherically as the force of capillarity controls the early growth (see 
Figure 3a and 3b). The volume of the bubble at 7 = 0.683 s and 1.967 s is 1.937 mm^ and 
5.473 mm3, respectively. By 7 = 2.927 s, the force of buoyancy becomes important on the 
bubble of volume 8.176 mm3, rjpg bubble translates upwards and begins to lose sphericity 
above the contact line (see Figure [^. 

As the bubble seeks to minimise its surface area at a given volume and due to the fact 
that the hydrostatic pressure varies linearly along the bubble, a neck begins to form in the 
free surface. This can be seen at 7 = 4.137 s (see Figure 3d) where the volume of the 
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Figure 4: The final computed solution of the case {rc,Oh,Q) = (0.106,2.24 x 10 ^,5.11 x 
10“®), as Tmin = rtoi, where i = 4.139 s and so td = 1.107 x 10®. 


bubble is now 11.518 mm®. Once the neck has formed, as there is a change in the sign 
of the longitudinal curvature at some point on the free surface, the pinch-off process can 
begin. The increasing difference in hydrostatic pressure between the apex and the base of 
the bubble results in an increase of capillary pressure at the neck which then drives the 
shrinking of the neck still further. 

Figure shows the final computed solution, as = rtoi, at i = 4.139 s. Here, the 
total volume of the bubble is 11.523 mm®, whilst the volume of the newly formed bubble, 
measured above the point of minimum neck radius is 11.502 mm®. Therefore, this final solu¬ 
tion corresponds to a dimensionless formation time of td = 1.107 x 10® and a dimensionless 
bubble volume of 14 = 0.5651. 

Compared to the overall generation of the bubble, which took 4.139 s, the pinch-off 
process takes place much faster in 1.68 x 10“^ s. Notably, while the experiments appear 
unable to capture the very small time scales in the final stages of pinch-off, this is no barrier 
for the computations. 

The evolution of the dimensional gas pressure during the simulation is shown by Figure]^ 
Recalling that the pressure field along the side boundary of the domain remains hydrostatic 
with p{R, 0, t) = 0, the gas pressure reaches its greatest value near the start of the simulation 
as the bubble resembles a hemisphere. The gas pressure then decreases by almost 80% until 
there is a very slight increase as t rtoi, which can only be seen under the appropriate 

magnification. 


14 






Figure 5: The evolution of the dimensional gas pressure Pg with dimensional simulation time 
i for the case of {rc,Oh,Q) = (0.106,2.24 x 10-3,5.11 x 10"®). 


5.2. Influence of Parameters 

Figure shows how the global characteristics of bubble formation, the dimensionless 
formation time (see Figure]^ and the dimensionless bubble volume (see Figure [hb|, 
depend upon the gas flow rate Q and the Ohnesorge number Oh for the dimensionless orihce 
radii of Tc = 0.1 and 1. These benchmark calculations fall into two regimes, the regimes of 
low gas flow rates and high gas flow rates, otherwise known as the ‘static’ and ‘dynamic’ 
regimes, respectively. These regimes are now considered separately. 


5.2.1. Regime of Low Gas Flow Rates 

For a relatively small gas flow rate, a variation in the Ohnesorge number has very little 
effect on the formation time (see Figure [6a)) and, consequently, very little effect on the bubble 
volume (see Figure [6b|. Therefore, for a given orihce radius, as Q —?• 0, then Vd —)■ 14, where 
Vc is the limiting bubble volume. Specihcally, 14 = 0.53 and 4.7 for Uc = 0.1 and 1.0, 
respectively. In other words, for a given orihce radius, simply decreasing the how rate does 
not lead to smaller bubble volumes in this regime. The force of buoyancy is simply not large 
enough to detach the bubble from the formation site until the bubble volume reaches 14. 

By assuming that the bubble remains spherical and balancing the forces of buoyancy and 
capillarity, Fritz (1935) (see Kumar and Kuloor ( 1970[ )) derived an expression for this limiting 
bubble volume which, when rescaled using the scales found in the problem formulation, is 
given by, 

Vf = 2'Krc. (14) 


Therefore Vp = 6.28 x IQ-^ and 6.28 for Tc = 0.1 and 1, respectively (see lines 1 and 2 
respectively on Figure [61^. This gives way to large relative errors of 18.5% and 33.6% for 
Tc = 0.1 and 1, respectively, when compared to the simulated limiting volume I 4 . 
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(a) A log-log plot of the dimensionless formation time td of the bubble versus the 
dimensionless flow rate Q for a variety of Ohnesorge numbers Oh. Line 1 is = 
0.53/(5 a-nd line 2 is = 4.7/(5. 



(b) A log-log plot of the dimensionless bubble volume Vd versus the dimensionless 
flow rate Q for a variety of Ohnesorge numbers Oh. Line 1 is Vp = 6.28 x 10“^, line 
2 is Vf = 6.28, line 3 is equation ( |17[ ) with Oh = 2.24 x 10“^ and line 4 is equation 

(@. 

Figure 6: A map of the parameter space where open symbols correspond to = 0.1 and 
hlled symbols correspond to Tc = 1. The various Ohnesorge numbers are represented by (o) 
Oh = 2.24 X 10-3, (>) Oh = lO-^, (□) Oh = lO'S (o) Oh = 1 and (<) Oh = 10. 
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(a) The case of {rc,Oh,Q) = (0.1,1,10 ^). (b) The case of {rc,Oh,Q) = (1,1,10“^). 

Figure 7: The formation of bubbles in the regime of low gas flow rates. The dynamic finite 
element simulations are given by the black lines whilst the circular symbols represent the 
hnite difference Young-Laplace solutions for quasi-static growth. 


Since the volume of the newly formed bubble is much greater than the volume of the 
residual bubble then td ^ VcjQ in this regime. Therefore, for Tc = 0.1 and 1, the respective 
formation times can be expressed as td = 0.53/Q and td = 4:.7/Q when Q < 10““^ (see lines 
1 and 2 respectively in Figure [^. 

As expected, the limiting bubble volume 14 increases with orifice radius. These con¬ 
clusions are in agreement with the literature where the regime of low gas flow rates is also 
known as the ‘static’ regime. As will now be shown below, much of the growth of the bubble 
can be described by a quasi-static approach that involves the Young-Laplace equation. 

Figure [^shows the simulated temporal evolution of a bubble (continuous black lines) for 
two cases in the regime of low gas flow rates, {rc,Oh,Q) = (0.1,1,10“^) and (1,1,10“®), 
respectively. In both cases, curve 1 is the initial solution and curve 4 shows the free surface 
when the neck first begins to develop. Curves 2 and 3 are equally spaced in time between 
curves 1 and 4, where it can be seen that, once again, the capillarity force dominates the 
early stages of bubble growth. This is particularly evident in the case of = 0.1 (see 
Figure 7a) where the bubble grows spherically for almost the entire formation time due to 
its smaller size than the case of = 1.0 (see Figure [7b|). As the pinch-off process continues, 
the bubble approaches break up and the hnal solution, as = rtoi, is given by curve 5 
where a bubble of volume 14 = 0.532 or 4.90 are formed for Vc = 0.1 or 1, respectively. 

Due to the relatively small liquid velocities associated with the initial bubble growth in 
these cases, the initial growth can be accurately described by a series of quasi-static profiles 
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that are found by solving the Young-Laplace equation and stipulating a small increase in 
volume from one solution to the next. 

Assuming that liquid velocities are negligible and rescaling, the dimensionless Young- 
Laplace equation is, 

Pg + Z = Ki + K2, (15) 


where Ki and K 2 are the respective dimensionless longitudinal and cross-sectional curvatures 
at each point on the free surface and the liquid pressure is assumed to be hydrostatic 
p(r, z, t) = —z. 

Fordham (|1948 ) and Thoroddsen et ah (2005) described a scheme that can be used to 


solve (15), which involves a system of first-order ordinary differential equations. The finite 
difference method can be used to solve this system and so for a given finite element profile, 
the corresponding quasi-static solution, with the same bubble volume as the simulation, 
can then be found. The uniform gas pressure in (15) is found as part of each quasi-static 
solution. 

These series of quasi-static profiles are also shown in Figure with the circular sym¬ 
bols. There is very good agreement between the hrst three quasi-static and hnite element 
solutions. However, as the bubble continues to grow, the quasi-static solutions become 
less accurate, as seen by the fourth solutions. Any further quasi-static solutions are wildly 
inaccurate. 

In summary, in the regime of low gas flow rates, the initial growth of the bubble is a 
quasi-static process that is governed by the gas, hydrostatic and capillary pressures. As 
the neck forms and pinch-off begins, bubble growth is an essentially dynamic process, and 
unsurprisingly, the quasi-static approach using the Young-Laplace equation (15) is inade¬ 
quate in describing the effects of the large velocities associated with the thinning of the neck 
together with the importance of inertia and viscosity in the liquid. 


5.2.2. Regime of High Gas Flow Rates 

Greater gas flow rates than those associated with the previous regime are now considered. 
For a given Ohnesorge number, as the flow rate is increased, it is now the formation time 
that tends to a limiting value tc (see Figure 6a). In other words, when the flow rate is 
sufficiently large, a bubble can not form quicker than this limiting time. This effect is most 
clear in the case of Tc = 1. As Q increases, the difference between the formation time of 
the simulations and the formation time associated with the regime of low gas flow rates 
{td = 0.53/Q and 4:.7/Q for = 0.1 and 1, respectively) increases. In agreement with the 
literature, this results in an increase in bubble volume Vd with flow rate (see Figure [hb|. 

For a given gas flow rate, the formation time increases with decreasing Ohnesorge number 
and therefore, due to the bubble inflating at a constant volumetric flow rate, the bubble 
volume increases with decreasing Ohnesorge number. This is once again more obvious in 
the case of Vc = 1.0 as the formation times tend to a limiting value. As the Ohnesorge 
number increases, the limiting formation time decreases. For Vc = 1, G ~ 1130, 260, 30, 7 
and 5 for Oh = 2.24 x 10“^, 10“^, 10“^, 1 and 10, respectively (see Figure 6a). The reason 
for this is the increased inertia in the liquid opposes the pinching of the neck which leads to 
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a prolonging of the formation time. There is very little difference in the results between an 
Ohnesorge number of 1 and 10, suggesting that inertial effects become negligible for Oh > 1. 
Once again, it is worth reiterating that the bubbles formed at a given flow rate from a larger 
orifice have a larger formation time and therefore a greater volume. 


Oguz and Prosperetti (1993) used a two-stage model which includes the forces of inertia. 


capillarity and buoyancy to derive a critical flow rate at which, in the case of an inviscid 
liquid, bubble formation transitions from the low flow rate regime to the high flow rate 
regime. By rescaling, this critical flow rate is given by, 

= (16/3)^/® TT Oh r®/®. 


For a flow rate Q < Qcr, Vd = Vp, whilst for Q > Qcr, it was found that, 

'_Q y/' 




OhJ 


, dvr / 9 

^ “ Y 


(16) 


(17) 


where the subscript 0 represents the inviscid regime. Oh —)■ 0. 

To compare ( p)^ with the results presented here, line 3 in Figure 6b is ([T^ with Oh = 
2.24 X 10“^, the smallest Ohnesorge number simulated. The simulated bubble volumes for 
Oh = 2.24 X 10“^ and Tc = 1.0 approach those values predicted by (17) as the flow rate 
increases. It is unclear whether the scaling law is valid for r,, = 0.1 because the greater flow 
rates required could not be computed by the numerical platform. 

Once again, the case of Oh = 2.24 x 10“^ is used to examine the critical flow rate given 
by (16), where Qcr = 1-37 x 10“^ and 9.3 x 10“^ for Tc = 0.1 and 1, respectively. These flow 
rates are also given by the points of intersection of line 3 by line 1 for Tc = 0.1 and line 2 


for Tc = 1 in Figure 6b It can be seen from the simulations that (16) overestimates the flow 


rate at which inviscid bubble formation transitions from the regime of low flow rates to the 
regime of high flow rates. 


In the case of a highly viscous liquid, Wong et ah (1998) used a spherical model and 


balanced the viscous and capillarity forces to derive an expression for the bubble volume 
where Tc —t 0. When rescaled. 


14, = A^Q^I\ As = (600ir/3) , 


(18) 


where the subscript cx) represents the highly viscous regime of Oh —>■ cx). The volumes 


predicted by (18) underestimate those given by the simulations for the most viscous case of 
Oh = 10 for both = 0.1 and 1 (see line 4 on Figure [6b|. 


It is possible to interpolate between the scaling laws and ( |I^ , for the respective 
limiting cases of Oh —)■ 0 and Oh —)■ cx, to £t a scaling law to the results for intermediate 
Ohnesorge numbers shown in Figure 6b The scaling law for bubble volume is assumed to 
take the form. 


Vk = 


Q6/5+C f{Oh) 

f{Oh) ■ 


f{Oh) = 


Oh^!^ 


Then, applying 14 —)■ Vq as Oh —)■ 0, and 14 —t 14c 


h Oh^O + 1 • 
as Oh —)■ cx, gives. 


a — 1/^1, b — A 2 IA\^ c — —9 ^ 42 / 20 , 
19 


(19) 



















Figure 8: Comparing the scaling law given by equation (19) for Vd in the high gas flow 
rate regime, depicted by the dashed lines, against the results presented in Figure ^ for 
Oh = 2.24 X 10-3 (line 1), Oh = IQ-^ (2), Oh = 10-\ (3), Oh = 1, (4) and Oh = 10 (5). 


where Vq and Ai are given by 0 and Voo and A 2 are given by 

For each Ohnesorge number investigated in this work. Figure [^shows that as the flow rate 
Q increases, the results from the simulations tend towards those predicted by the scaling law 
(19) as Q —)■ cxD, but do not give an accurate representation for moderate Q. Therefore, whilst 
scaling laws may give a valid approximation of the global characteristics of bubble formation 
in certain regimes, simulations are required to obtain a more accurate representation of the 
entire parameter space. 

Figures and IT show the evolution of the bubble for a cross-section of the parameter 
space in the regime of high gas flow rates for Vc = 0.1 and 1.0, respectively. As discussed 
previously, the cases of Oh = 1 and Oh = 10 are very similar and so the case of Oh = 10 is 
not given in Figures and 

The plots in Figure |^each show six curves. Curve 1 is the initial solution, curve 5 is 
where the bubble hrst obtains a neck and curve 6 is the hnal solution as Curves 

2, 3 and 4 are equally spaced in time between curves 1 and 5. The plots in Figure IT have 
an extra curve. Here, curve 4 is where the neck begins to develop and curve 7 is the hnal 
solution, curves 2 and 3 are equally spaced in time between curves 1 and 4, and curves 5 
and 6 show the bubble when = 2/3 Tc and = 1/3 Tc, respectively. 

Due to the greater velocity imparted in the liquid from the larger how rates and because 
inertia cannot be ignored at these how rates, the aforementioned quasi-static approach, 
using the Young-Laplace equation (15), to compute free surface prohles, is valid for a much 
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(a) Oh = 2.24 x IQ-^, (b) Oh = IQ-^, 

Q = 10-3. Q = 10-3. 


(c) Oh = 10 
Q = 10-3. 


(d) Oh = 1, 
Q = 10-3. 



(e) Oh = 2.24 X 10-3, 
g = 5 X 10-1 



(f) Oh = 10-2, 
g = 5 X 10-1 


z 

1.5 


0.5 



(g) Oh = 10 3, 
g = 5 X 10-1 


(h) Oh = 1, 
g = 5 X 10-1 






(i) Oh = 2.24 X 10-3, (j) Oh = 10-^, 

g = 10-1 g = 10-1 


(k) Oh = 10-3, 

g = 10-1 


(1) Oh = 1, 

g = 10-1 


Figure 9: Temporal evolution of the free surface for various Ohnesorge numbers Oh and flow 
rates Q for an orifice of radius = 0.1. 
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(a) Oh = 2.24x 10“3, 
Q = 10-2. 



(e) Oh = 2.24 X 10-3, 
g = 5 X 10-3. 



(i) Oh = 2.24 X 10-3, 

g = 10-3. 



(b) Oh = 10-2, 

g = 10 - 2 . 


(c) Oh = 10 

g = 10 - 2 . 


(d) Oh = 1, 

g = 10 - 2 . 



(f) Oh = 10-2, 
g = 5 X 10-3. 


(g) Oh = 10 3, 
g = 5 X 10-3. 


(h) Oh = 1, 
g = 5 X 10-3. 



(j) Oh = 10 2, 

g = 10-3. 


(k) Oh = 10-3, 

g = 10-3. 


(1) Oh = 1, 

g = 10-3. 


Figure 10: Temporal evolution of the free surface for various Ohnesorge numbers Oh and 
flow rates Q for an orifice of radius Tc = 1.0. 
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smaller proportion of the overall formation process than was seen in Figure for bubble 
formation under low gas flow rates. For greater flow rates, the bubble formation process is 
almost entirely dynamic and therefore this regime of high gas flow rates is also known as 
the ‘dynamic’ regime. 

A major difference between bubbles forming from orifices of different radii is the increased 
sphericity of those bubbles for = 0.1 due to the larger surface-to-volume ratio where 
surface effects are more important for smaller systems. Also the bubbles generated in the 
case of Tc = 0.1 are much larger relative to the orihce radius than those in the case of 
Tc = 1.0, which is to be expected since Vc = 1.0 corresponds to the dimensional case where 
the orihce radius is equal to the length scale. 

Let td = tg + tn, where tg is the time taken by the initial growth stage of the bubble, 
until the longitudinal curvature of the free surface changes sign at a point as the neck forms, 
and tn is the pinch-off time, taken to be the time taken for = f'toi once the neck has 
formed. 

For a given how rate, the ratio tn/td increases with decreasing Ohnesorge number. The 
inhuence of inertia in prolonging the formation time by opposing the thinning of the neck 
of the bubble is illustrated by Figures [M| and for Vc = 0.1, and by Figures 1^ and 1^ for 
Tc = 1.0, where the times taken for the necks to form are almost equal, yet td increases with 
decreasing Ohnesorge number. This is further highlighted by the rapid vertical displacement 
of the apex of the bubble seen in Figure lOd , where, in order to keep the how rate constant 
a.t Q = 10“^, the apex must rise quickly to make up the volume lost due to this rapid 
pinch-oh. 

For a given Ohnesorge number, as the how rate increases, the majority of the increase in 
volume is accounted for by the increased width of the bubbles, see for example. Figures |9l] 
and 9a for Vc = 0.1, and Figures lOi and 10a, for Vc = 1. The ratio tn/td initially increases 
with increasing how rate as seen in Figure]^ and IT As the how rate continues to increase 
tn/td reaches a maximum before decreasing with increasing how rate. 

Figure 11 shows the inhuence of the Ohnesorge number on the pinch-oh process. The 
outermost and innermost curves on each plot show the free surface as = 'f'c/‘l and 
'f'min = Ttoh respectively. The eight curves in between these are equally spaced in time with 
time step At = 17.6, 0.67 and 0.25 for Oh = 2.24 x 10“^, 10“^ and 10, respectively. 

The spacing of the curves show that the pinch-oh of the bubble accelerates as Tmin —t 0, 
whilst the radius of curvature at the point on the free surface which represents the minimum 
neck radius increases with increasing Ohnesorge number. 

Finally, for larger Ohnesorge numbers at sufficiently high how rates when Vc = 1, the 
liquid phase enters the mouth of the orifice towards the end of the formation process (see 
Figure [I^. The gas phase with spatially uniform pressure Pg is present not just above the 
orifice but below it and the solid surface too. Therefore there is no restriction on the free 
surface entering the mouth of the orifice if the normal stress boundary condition dictates 
this. 
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(a) Oh = 2.24 x 10 ^ and (b) Oh = 10 ^ and At = (c) Oh = 10 and At = 

At = 17.6. 0.67. 0.25. 

Figure 11: The influence of Ohnesorge number on the temporal evolution of the neck of 
the bubble during the pinch-off process for the case Tc = 1 and Q = 10“^. The outermost 
and innermost curves represent the respective cases of = ''"c/2 and = Ooi- The 
intermediate curves are equally spaced in time, with time step At. 


6. Concluding Remarks 

It was shown that the numerical platform based on the finite element method produced 
benchmark calculations for the formation of a single bubble from a submerged orifice. The 
global characteristics of the process, the formation time and bubble volume, were found for 
a wide range of orifice radii, Ohnesorge numbers and volumetric gas flow rates. 

The results were seen to agree very well with experiments, as the finite element method 
allows the free surface of the bubble to be tracked explicitly and accurately, even when the 
liquid phase enters the mouth of the orifice, while the computational mesh can be refined 
in the regions of the flow domain that require extra precision. The simulations could also 
capture the very small time scales associated with the pinch-off of the bubble whilst the 
experiments were unable to do so. 

In agreement with the literature, two regimes of bubble formation were identified, the 
low gas low rate ‘static’ regime and the high gas flow rate ‘dynamic’ regime. The benchmark 
calculations were used to validate a range of scaling laws proposed in the literature. 

The liquid phase was assumed to perfectly wet the solid surface throughout the entire 
formation process, and so the three phase solid-liquid-gas contact line remained pinned 
to the rim of the orifice, whilst the gas was assumed to be inviscid. The model could 
be extended by relaxing these assumptions in order to investigate the influence of partial 
wetting conditions, whereby the contact line is free to move along the solid surface, and the 
viscosity of the gas on the bubble formation process. 

The process of the break up of the bubble, where the topology of the domain occupied 
by the gas changes, presents a major problem for the modelling. In the present work, it is 
assumed to have a local effect which is an acceptable approximation for the relatively large 
bubbles considered here. In a general case, however, especially in microfluidics, it becomes 
necessary to model the break up process in a singularity-free way. Currently, this problem 
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Figure 12: Evidence of the liquid phase entering the mouth of the orihce as = rtoi for 
the case (rc, Oh, Q) = (1,10, 7.5) where td = 6.596 and Vd = 50.12. 


is outstanding and requires further research. 


Appendix A. Rescaling the Resnlts 

In order to convert the results into the dimensionless system used in some works on 
bubble dynamics, recall that fc = Lvc and Q = L^UQ are the respective dimensional orifice 
radius and volumetric gas flow rate. In some works fc and Q/f1 are used as the length and 
velocity scales rather than those used here of L = (J jpg and U = a/p, respectively. 

Rescaling in this manner leads to a different group of dimensionless parameters, which 
will be denoted with hats. Rather than (rc. Oh, Q), there is [Bo, Oh, hFej or {^Bo, Oh, Ca^ , 

where Bo = pgfl/a. Oh = p,/y/pfc<J, We = pQ'^/af^ and Ca = pQ/fla, are the rescaled 
Bond, Ohnesorge, Weber and capillary nnmbers. The rescaled dimensionless formation time 
and volume are denoted by td and Vd- Then, to convert the parameters used here into the 
rescaled parameters, we have. 
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Bo 
Oh 
We 
Ca 
td 
Vd 

and the inverse is, 

Oh 

Q 

td 

Vd 
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